multinichenetr 2.0.0
In this vignette, you can learn how to perform a MultiNicheNet analysis to compare cell-cell communication between conditions of interest. A MultiNicheNet analysis can be performed if you have multi-sample, multi-condition/group single-cell data. We strongly recommend having at least 4 samples in each of the groups/conditions you want to compare. With less samples, the benefits of performing a pseudobulk-based DE analysis are less clear.
As input you need a SingleCellExperiment object containing at least the raw count matrix and metadata providing the following information for each cell: the group, sample and cell type.
As example expression data of interacting cells for this vignette, we will here use scRNAseq data from breast cancer biopsies of patients receiving anti-PD1 immune-checkpoint blockade therapy. Bassez et al. collected from each patient one tumor biopsy before anti-PD1 therapy (“pre-treatment”) and one during subsequent surgery (“on-treatment”) A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Based on additional scTCR-seq results, they identified one group of patients with clonotype expansion as response to the therapy (“E”) and one group with only limited or no clonotype expansion (“NE”).
We will use MultiNicheNet to explore immune cell crosstalk specfic for Expander patients compared to Non-expander patients, focusing on the data BEFORE treatment.
In this vignette, we will first prepare the MultiNicheNet core analysis, then run the several steps in the MultiNicheNet core analysis, and finally interpret the output.
library(SingleCellExperiment)
library(nichenetr)
library(multinichenetr)
library(tidyverse)
MultiNicheNet builds upon the NicheNet framework and uses the same prior knowledge networks (ligand-receptor network and ligand-target matrix, currently v2 version).
The Nichenet v2 networks and matrices for both mouse and human can be downloaded from Zenodo .
We will read these object in for human because our expression data is of human patients.
Gene names are here made syntactically valid via make.names() to avoid the loss of genes (eg H2-M3) in downstream visualizations.
organism = "human"
options(timeout = 120)
if(organism == "human"){
lr_network_all =
readRDS(url(
"https://zenodo.org/record/10229222/files/lr_network_human_allInfo_30112033.rds"
)) %>%
mutate(
ligand = convert_alias_to_symbols(ligand, organism = organism),
receptor = convert_alias_to_symbols(receptor, organism = organism))
lr_network_all = lr_network_all %>%
mutate(ligand = make.names(ligand), receptor = make.names(receptor))
lr_network = lr_network_all %>%
distinct(ligand, receptor)
ligand_target_matrix = readRDS(url(
"https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds"
))
colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>%
convert_alias_to_symbols(organism = organism) %>% make.names()
rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>%
convert_alias_to_symbols(organism = organism) %>% make.names()
lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix))
ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()]
} else if(organism == "mouse"){
lr_network_all = readRDS(url(
"https://zenodo.org/record/10229222/files/lr_network_mouse_allInfo_30112033.rds"
)) %>%
mutate(
ligand = convert_alias_to_symbols(ligand, organism = organism),
receptor = convert_alias_to_symbols(receptor, organism = organism))
lr_network_all = lr_network_all %>%
mutate(ligand = make.names(ligand), receptor = make.names(receptor))
lr_network = lr_network_all %>%
distinct(ligand, receptor)
ligand_target_matrix = readRDS(url(
"https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds"
))
colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>%
convert_alias_to_symbols(organism = organism) %>% make.names()
rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>%
convert_alias_to_symbols(organism = organism) %>% make.names()
lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix))
ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()]
}
In this vignette, we will load in a subset of the scRNAseq data of the BRCA . For the sake of demonstration, this subset only contains 3 cell types. These celltypes are some of the cell types that were found to be most interesting related to BRCA according to Hoste et al.
If you start from a Seurat object, you can convert it easily to a SingleCellExperiment object via sce = Seurat::as.SingleCellExperiment(seurat_obj, assay = "RNA").
Because the NicheNet 2.0. networks are in the most recent version of the official gene symbols, we will make sure that the gene symbols used in the expression data are also updated (= converted from their “aliases” to official gene symbols). Afterwards, we will make them again syntactically valid.
options(timeout = 500)
sce = readRDS(url(
"https://zenodo.org/record/8010790/files/sce_subset_breastcancer.rds"
))
sce = alias_to_symbol_SCE(sce, "human") %>% makenames_SCE() # why?: NicheNet-v2 in recent gene symbols
In this step, we will formalize our research question into MultiNicheNet input arguments.
In this case study, we want to study differences in cell-cell communication patterns between “expander” BRCA patients (PreE) before therapy and “non-expander” patients before therapy (PreNE). The meta data columns that indicate this status (=group/condition of interest) is expansion_timepoint.
Cell type annotations are indicated in the subType column, and the sample is indicated by the sample_id column.
If your cells are annotated in multiple hierarchical levels, we recommend using a relatively high level in the hierarchy. This for 2 reasons: 1) MultiNicheNet focuses on differential expression and not differential abundance, and 2) there should be sufficient cells per sample-celltype combination (see later).
sample_id = "sample_id"
group_id = "expansion_timepoint"
celltype_id = "subType"
Important: It is required that each sample-id is uniquely assigned to only one condition/group of interest. See the vignettes about paired and multifactorial analysis to see how to define your analysis input when you have multiple samples (and conditions) per patient.
If you would have batch effects or covariates you can correct for, you can define this here as well. However, this is not applicable to this dataset. Therefore we will use the following NA settings:
covariates = NA
batches = NA
Important: for categorical covariates and batches, there should be at least one sample for every group-batch combination. If one of your groups/conditions lacks a certain level of your batch, you won’t be able to correct for the batch effect because the model is then not able to distinguish batch froPreE group/condition effects.
Important: The column names of group, sample, cell type, batches and covariates should be syntactically valid (run make.names)
Important: All group, sample, cell type, batch and covariate names should be syntactically valid as well (run make.names) (eg through SummarizedExperiment::colData(sce)$ShortID = SummarizedExperiment::colData(sce)$ShortID %>% make.names())
Here, we want to know which cell-cell communication patterns are specific for the PreE group versus the PreNE group and vice versa.
To perform this comparison, we need to set the following contrasts:
contrasts_oi = c("'PreE-PreNE','PreNE-PreE'")
Very Important Note the format to indicate the contrasts! This formatting should be adhered to very strictly, and white spaces are not allowed! Check ?get_DE_info for explanation about how to define this well. The most important points are that:
each contrast is surrounded by single quotation marks
contrasts are separated by a comma without any white space
*all contrasts together are surrounded by double quotation marks.
If you compare against two groups, you should divide by 2 (as demonstrated in another vignette), if you compare against three groups, you should divide by 3 and so on.
For downstream visualizations and linking contrasts to their main condition, we also need to run the following: This is necessary because we will also calculate cell-type+condition specificity of ligands and receptors.
contrast_tbl = tibble(contrast =
c("PreE-PreNE","PreNE-PreE"),
group = c("PreE","PreNE"))
If you want to focus the analysis on specific cell types (e.g. because you know which cell types reside in the same microenvironments based on spatial data), you can define this here. If you have sufficient computational resources and no specific idea of cell-type colocalzations, we recommend to consider all cell types as potential senders and receivers. Later on during analysis of the output it is still possible to zoom in on the cell types that interest you most, but your analysis is not biased to them.
Here we will consider all cell types in the data:
senders_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()
receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()
sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in%
c(senders_oi, receivers_oi)
]
In case you would have samples in your data that do not belong to one of the groups/conditions of interest, we recommend removing them and only keeping conditions of interst. Why? To determine expressed genes in the gene filtering step, it is best to only keep conditions that are of direct interest.
conditions_keep = c("PreE", "PreNE")
sce = sce[, SummarizedExperiment::colData(sce)[,group_id] %in%
conditions_keep
]
Now we will run the core of a MultiNicheNet analysis. This analysis consists of the following steps:
Following these steps, one can optionally * 7. Calculate the across-samples expression correlation between ligand-receptor pairs and target genes * 8. Prioritize communication patterns involving condition-specific cell types through an alternative prioritization scheme
After these steps, the output can be further explored as we will demonstrate in the “Downstream analysis of the MultiNicheNet output” section.
In this step we will calculate and visualize cell type abundances. This will give an indication about which cell types will be retained in the analysis, and which cell types will be filtered out.
Since MultiNicheNet will infer group differences at the sample level for each cell type (currently via Muscat - pseudobulking + EdgeR), we need to have sufficient cells per sample of a cell type, and this for all groups. In the following analysis we will set this minimum number of cells per cell type per sample at 10. Samples that have less than min_cells cells will be excluded from the analysis for that specific cell type.
min_cells = 10
We recommend using min_cells = 10, except for datasets with several lowly abundant cell types of interest. For those datasets, we recommend using min_cells = 5.
abundance_info = get_abundance_info(
sce = sce,
sample_id = sample_id, group_id = group_id, celltype_id = celltype_id,
min_cells = min_cells,
senders_oi = senders_oi, receivers_oi = receivers_oi,
batches = batches
)
First, we will check the cell type abundance diagnostic plots.
The first plot visualizes the number of cells per celltype-sample combination, and indicates which combinations are removed during the DE analysis because there are less than min_cells in the celltype-sample combination.
abundance_info$abund_plot_sample
The red dotted line indicates the required minimum of cells as defined above in min_cells. We can see here that some sample-celltype combinations are left out. For the DE analysis in the next step, only cell types will be considered if there are at least two samples per group with a sufficient number of cells. But as we can see here: all cell types will be considered for the analysis and there are no condition-specific cell types.
Important: Based on the cell type abundance diagnostics, we recommend users to change their analysis settings if required (eg changing cell type annotation level, batches, …), before proceeding with the rest of the analysis. If too many celltype-sample combinations don’t pass this threshold, we recommend to define your cell types in a more general way (use one level higher of the cell type ontology hierarchy) (eg TH17 CD4T cells –> CD4T cells) or use min_cells = 5 if this would not be possible.
Running the following block of code can help you determine which cell types are condition-specific and which cell types are absent.
sample_group_celltype_df = abundance_info$abundance_data %>%
filter(n > min_cells) %>%
ungroup() %>%
distinct(sample_id, group_id) %>%
cross_join(
abundance_info$abundance_data %>%
ungroup() %>%
distinct(celltype_id)
) %>%
arrange(sample_id)
abundance_df = sample_group_celltype_df %>% left_join(
abundance_info$abundance_data %>% ungroup()
)
abundance_df$n[is.na(abundance_df$n)] = 0
abundance_df$keep[is.na(abundance_df$keep)] = FALSE
abundance_df_summarized = abundance_df %>%
mutate(keep = as.logical(keep)) %>%
group_by(group_id, celltype_id) %>%
summarise(samples_present = sum((keep)))
celltypes_absent_one_condition = abundance_df_summarized %>%
filter(samples_present == 0) %>% pull(celltype_id) %>% unique()
# find truly condition-specific cell types by searching for cell types
# truely absent in at least one condition
celltypes_present_one_condition = abundance_df_summarized %>%
filter(samples_present >= 2) %>% pull(celltype_id) %>% unique()
# require presence in at least 2 samples of one group so
# it is really present in at least one condition
condition_specific_celltypes = intersect(
celltypes_absent_one_condition,
celltypes_present_one_condition)
total_nr_conditions = SummarizedExperiment::colData(sce)[,group_id] %>%
unique() %>% length()
absent_celltypes = abundance_df_summarized %>%
filter(samples_present < 2) %>%
group_by(celltype_id) %>%
count() %>%
filter(n == total_nr_conditions) %>%
pull(celltype_id)
print("condition-specific celltypes:")
## [1] "condition-specific celltypes:"
print(condition_specific_celltypes)
## character(0)
print("absent celltypes:")
## [1] "absent celltypes:"
print(absent_celltypes)
## character(0)
Absent cell types will be filtered out, condition-specific cell types can be filtered out if you as a user do not want to run the alternative workflow for condition-specific cell types in the optional step 8 of the core MultiNicheNet analysis.
analyse_condition_specific_celltypes = FALSE
if(analyse_condition_specific_celltypes == TRUE){
senders_oi = senders_oi %>% setdiff(absent_celltypes)
receivers_oi = receivers_oi %>% setdiff(absent_celltypes)
} else {
senders_oi = senders_oi %>%
setdiff(union(absent_celltypes, condition_specific_celltypes))
receivers_oi = receivers_oi %>%
setdiff(union(absent_celltypes, condition_specific_celltypes))
}
sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in%
c(senders_oi, receivers_oi)
]
Before running the DE analysis, we will determine which genes are not sufficiently expressed and should be filtered out.
We will perform gene filtering based on a similar procedure as used in edgeR::filterByExpr. However, we adapted this procedure to be more interpretable for single-cell datasets.
For each cell type, we will consider genes expressed if they are expressed in at least a min_sample_prop fraction of samples in the condition with the lowest number of samples. By default, we set min_sample_prop = 0.50, which means that genes should be expressed in at least 4 samples if the group with lowest nr. of samples has 9 samples like this dataset.
min_sample_prop = 0.50
But how do we define which genes are expressed in a sample? For this we will consider genes as expressed if they have non-zero expression values in a fraction_cutoff fraction of cells of that cell type in that sample. By default, we set fraction_cutoff = 0.05, which means that genes should show non-zero expression values in at least 5% of cells in a sample.
fraction_cutoff = 0.05
We recommend using these default values unless there is specific interest in prioritizing (very) weakly expressed interactions. In that case, you could lower the value of fraction_cutoff. We explicitly recommend against using fraction_cutoff > 0.10.
Now we will calculate the information required for gene filtering with the following command:
frq_list = get_frac_exprs(
sce = sce,
sample_id = sample_id, celltype_id = celltype_id, group_id = group_id,
batches = batches,
min_cells = min_cells,
fraction_cutoff = fraction_cutoff, min_sample_prop = min_sample_prop)
## [1] "Samples are considered if they have more than 10 cells of the cell type of interest"
## [1] "Genes with non-zero counts in at least 5% of cells of a cell type of interest in a particular sample will be considered as expressed in that sample."
## [1] "Genes expressed in at least 4.5 samples will considered as expressed in the cell type: CD4T"
## [1] "Genes expressed in at least 4.5 samples will considered as expressed in the cell type: Fibroblast"
## [1] "Genes expressed in at least 4.5 samples will considered as expressed in the cell type: macrophages"
## [1] "7945 genes are considered as expressed in the cell type: CD4T"
## [1] "9524 genes are considered as expressed in the cell type: Fibroblast"
## [1] "9138 genes are considered as expressed in the cell type: macrophages"
Now only keep genes that are expressed by at least one cell type:
genes_oi = frq_list$expressed_df %>%
filter(expressed == TRUE) %>% pull(gene) %>% unique()
sce = sce[genes_oi, ]
After filtering out absent cell types and genes, we will continue the analysis by calculating the different prioritization criteria that we will use to prioritize cell-cell communication patterns.
First, we will determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type. The function process_abundance_expression_info will link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types. This will later on allow us to define the cell-type specicificy criteria for ligands and receptors.
abundance_expression_info = process_abundance_expression_info(
sce = sce,
sample_id = sample_id, group_id = group_id, celltype_id = celltype_id,
min_cells = min_cells,
senders_oi = senders_oi, receivers_oi = receivers_oi,
lr_network = lr_network,
batches = batches,
frq_list = frq_list,
abundance_info = abundance_info)
Normalized pseudobulk expression values per gene/celltype/sample can be inspected by:
abundance_expression_info$celltype_info$pb_df %>% head()
## # A tibble: 6 × 4
## gene sample pb_sample celltype
## <chr> <chr> <dbl> <fct>
## 1 A1BG BIOKEY_10Pre 5.09 CD4T
## 2 A1BG.AS1 BIOKEY_10Pre 3.64 CD4T
## 3 A2M BIOKEY_10Pre 2.92 CD4T
## 4 A4GALT BIOKEY_10Pre 0.526 CD4T
## 5 AAAS BIOKEY_10Pre 4.73 CD4T
## 6 AADACL2.AS1 BIOKEY_10Pre 0 CD4T
An average of these sample-level expression values per condition/group can be inspected by:
abundance_expression_info$celltype_info$pb_df_group %>% head()
## # A tibble: 6 × 4
## # Groups: group, celltype [1]
## group celltype gene pb_group
## <chr> <chr> <chr> <dbl>
## 1 PreE CD4T A1BG 5.39
## 2 PreE CD4T A1BG.AS1 3.35
## 3 PreE CD4T A2M 2.25
## 4 PreE CD4T A4GALT 0.553
## 5 PreE CD4T AAAS 4.31
## 6 PreE CD4T AADACL2.AS1 0
Inspecting these values for ligand-receptor interactions can be done by:
abundance_expression_info$sender_receiver_info$pb_df %>% head()
## # A tibble: 6 × 8
## sample sender receiver ligand receptor pb_ligand pb_receptor
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 BIOKEY_6Pre macrophages macrophages HLA.DMA CD74 10.9 15.6
## 2 BIOKEY_31Pre macrophages macrophages MIF CD74 12.1 13.9
## 3 BIOKEY_24Pre macrophages macrophages HLA.DMA CD74 10.8 15.3
## 4 BIOKEY_14Pre macrophages macrophages HLA.DMA CD74 10.8 15.3
## 5 BIOKEY_4Pre macrophages macrophages HLA.DMA CD74 10.8 15.3
## 6 BIOKEY_27Pre macrophages macrophages HLA.DMA CD74 10.9 15.1
## # ℹ 1 more variable: ligand_receptor_pb_prod <dbl>
abundance_expression_info$sender_receiver_info$pb_df_group %>% head()
## # A tibble: 6 × 8
## # Groups: group, sender [5]
## group sender receiver ligand receptor pb_ligand_group pb_receptor_group
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 PreE macrophages macrophag… MIF CD74 10.1 14.5
## 2 PreNE macrophages macrophag… HLA.D… CD74 10.0 14.6
## 3 PreNE Fibroblast macrophag… MIF CD74 10.0 14.6
## 4 PreE macrophages macrophag… HLA.D… CD74 9.87 14.5
## 5 PreE Fibroblast macrophag… MIF CD74 9.72 14.5
## 6 PreE CD4T macrophag… MIF CD74 9.70 14.5
## # ℹ 1 more variable: ligand_receptor_pb_prod_group <dbl>
In this step, we will perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest (as formalized by the contrasts_oi). Based on this analysis, we later can define the levels of differential expression of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver cell types (which will be used for the ligand activity analysis).
We will apply pseudobulking followed by EdgeR to perform multi-condition multi-sample differential expression (DE) analysis (also called ‘differential state’ analysis by the developers of Muscat).
DE_info = get_DE_info(
sce = sce,
sample_id = sample_id, group_id = group_id, celltype_id = celltype_id,
batches = batches, covariates = covariates,
contrasts_oi = contrasts_oi,
min_cells = min_cells,
expressed_df = frq_list$expressed_df)
## [1] "DE analysis is done:"
## [1] "included cell types are:"
## [1] "CD4T" "Fibroblast" "macrophages"
Check DE output information in table with logFC and p-values for each gene-celltype-contrast:
DE_info$celltype_de$de_output_tidy %>% head()
## # A tibble: 6 × 9
## gene cluster_id logFC logCPM F p_val p_adj.loc p_adj contrast
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 A1BG CD4T -0.343 5.69 3.54 0.0677 0.442 0.442 PreE-PreNE
## 2 AAAS CD4T 0.00365 4.43 0.000335 0.985 1 1 PreE-PreNE
## 3 AAGAB CD4T 0.024 5.08 0.02 0.888 1 1 PreE-PreNE
## 4 AAK1 CD4T -0.145 7.11 0.726 0.4 0.839 0.839 PreE-PreNE
## 5 AAMDC CD4T -0.282 3.8 1.02 0.319 0.786 0.786 PreE-PreNE
## 6 AAMP CD4T 0.0023 5.96 0.000186 0.989 1 1 PreE-PreNE
Evaluate the distributions of p-values:
DE_info$hist_pvals
These distributions look fine (uniform distribution, except peak at p-value <= 0.05), so we will continue using these regular p-values. In case these p-value distributions look irregular, you can estimate empirical p-values as we will demonstrate in another vignette.
empirical_pval = FALSE
if(empirical_pval == TRUE){
DE_info_emp = get_empirical_pvals(DE_info$celltype_de$de_output_tidy)
celltype_de = DE_info_emp$de_output_tidy_emp %>% select(-p_val, -p_adj) %>%
rename(p_val = p_emp, p_adj = p_adj_emp)
} else {
celltype_de = DE_info$celltype_de$de_output_tidy
}
To end this step, we will combine the DE information of senders and receivers by linking their ligands and receptors together based on the prior knowledge ligand-receptor network.
sender_receiver_de = combine_sender_receiver_de(
sender_de = celltype_de,
receiver_de = celltype_de,
senders_oi = senders_oi,
receivers_oi = receivers_oi,
lr_network = lr_network
)
sender_receiver_de %>% head(20)
## # A tibble: 20 × 12
## contrast sender receiver ligand receptor lfc_ligand lfc_receptor
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 PreE-PreNE CD4T CD4T CCL4L2 CCR1 4.53 2.11
## 2 PreE-PreNE macrophages CD4T TGM2 ADGRG1 2 3.7
## 3 PreE-PreNE Fibroblast CD4T MMP1 ITGA2 3.41 1.78
## 4 PreE-PreNE CD4T macrophages GZMB IGF2R 4.75 0.328
## 5 PreE-PreNE CD4T macrophages CCL4L2 CCR5 4.53 0.508
## 6 PreE-PreNE CD4T Fibroblast GZMB MCL1 4.75 0.279
## 7 PreE-PreNE CD4T CD4T CCL4L2 CCR5 4.53 0.39
## 8 PreE-PreNE macrophages CD4T CCL13 CCR1 2.79 2.11
## 9 PreE-PreNE CD4T macrophages CSF2 CSF2RB 4.1 0.737
## 10 PreE-PreNE CD4T macrophages CCL4L2 CCR1 4.53 0.266
## 11 PreE-PreNE CD4T Fibroblast GZMB IGF2R 4.75 0.0425
## 12 PreE-PreNE CD4T CD4T GZMB MCL1 4.75 0.0195
## 13 PreE-PreNE Fibroblast CD4T TGM2 ADGRG1 1.04 3.7
## 14 PreE-PreNE macrophages CD4T CCL5 CCR1 2.61 2.11
## 15 PreE-PreNE CD4T CD4T GZMB IGF2R 4.75 -0.0636
## 16 PreE-PreNE CD4T macrophages CSF2 IL3RA 4.1 0.49
## 17 PreE-PreNE CD4T macrophages GZMB MCL1 4.75 -0.192
## 18 PreNE-PreE CD4T CD4T CXCL14 CXCR4 3.94 0.537
## 19 PreNE-PreE CD4T macrophages CXCL14 CXCR4 3.94 0.474
## 20 PreE-PreNE macrophages CD4T CCL5 SDC4 2.61 1.61
## # ℹ 5 more variables: ligand_receptor_lfc_avg <dbl>, p_val_ligand <dbl>,
## # p_adj_ligand <dbl>, p_val_receptor <dbl>, p_adj_receptor <dbl>
In this step, we will predict NicheNet ligand activities and NicheNet ligand-target links based on these differential expression results. We do this to prioritize interactions based on their predicted effect on a receiver cell type. We will assume that the most important group-specific interactions are those that lead to group-specific gene expression changes in a receiver cell type.
Similarly to base NicheNet (https://github.com/saeyslab/nichenetr), we use the DE output to create a “geneset of interest”: here we assume that DE genes within a cell type may be DE because of differential cell-cell communication processes. In the ligand activity prediction, we will assess the enrichment of target genes of ligands within this geneset of interest. In case high-probabiliy target genes of a ligand are enriched in this set compared to the background of expressed genes, we predict that this ligand may have a high activity.
Because the ligand activity analysis is an enrichment procedure, it is important that this geneset of interest should contain a sufficient but not too large number of genes. The ratio geneset_oi/background should ideally be between 1/200 and 1/10 (or close to these ratios).
To determine the genesets of interest based on DE output, we need to define some logFC and/or p-value thresholds per cell type/contrast combination. In general, we recommend inspecting the nr. of DE genes for all cell types based on the default thresholds and adapting accordingly. By default, we will apply the p-value cutoff on the normal p-values, and not on the p-values corrected for multiple testing. This choice was made because most multi-sample single-cell transcriptomics datasets have just a few samples per group and we might have a lack of statistical power due to pseudobulking. But, if the smallest group >= 20 samples, we typically recommend using p_val_adj = TRUE. When the biological difference between the conditions is very large, we typically recommend increasing the logFC_threshold and/or using p_val_adj = TRUE.
We will first inspect the geneset_oi-vs-background ratios for the default tresholds:
logFC_threshold = 0.50
p_val_threshold = 0.05
p_val_adj = FALSE
geneset_assessment = contrast_tbl$contrast %>%
lapply(
process_geneset_data,
celltype_de, logFC_threshold, p_val_adj, p_val_threshold
) %>%
bind_rows()
geneset_assessment
## # A tibble: 6 × 12
## cluster_id n_background n_geneset_up n_geneset_down prop_geneset_up
## <chr> <int> <int> <int> <dbl>
## 1 CD4T 7945 409 230 0.0515
## 2 Fibroblast 9524 219 95 0.0230
## 3 macrophages 9138 576 597 0.0630
## 4 CD4T 7945 230 409 0.0289
## 5 Fibroblast 9524 95 219 0.00997
## 6 macrophages 9138 597 576 0.0653
## # ℹ 7 more variables: prop_geneset_down <dbl>, in_range_up <lgl>,
## # in_range_down <lgl>, contrast <chr>, logFC_threshold <dbl>,
## # p_val_threshold <dbl>, adjusted <lgl>
We can see here that for all cell type / contrast combinations, all geneset/background ratio’s are within the recommended range (in_range_up and in_range_down columns). When these geneset/background ratio’s would not be within the recommended ranges, we should interpret ligand activity results for these cell types with more caution, or use different thresholds (for these or all cell types).
For the sake of demonstration, we will also calculate these ratio’s in case we would use the adjusted p-value as threshold.
geneset_assessment_adjustedPval = contrast_tbl$contrast %>%
lapply(
process_geneset_data,
celltype_de, logFC_threshold, p_val_adj = TRUE, p_val_threshold
) %>%
bind_rows()
geneset_assessment_adjustedPval
## # A tibble: 6 × 12
## cluster_id n_background n_geneset_up n_geneset_down prop_geneset_up
## <chr> <int> <int> <int> <dbl>
## 1 CD4T 7945 156 36 0.0196
## 2 Fibroblast 9524 6 0 0.000630
## 3 macrophages 9138 40 24 0.00438
## 4 CD4T 7945 36 156 0.00453
## 5 Fibroblast 9524 0 6 0
## 6 macrophages 9138 24 40 0.00263
## # ℹ 7 more variables: prop_geneset_down <dbl>, in_range_up <lgl>,
## # in_range_down <lgl>, contrast <chr>, logFC_threshold <dbl>,
## # p_val_threshold <dbl>, adjusted <lgl>
We can see here that for most cell type / contrast combinations, the geneset/background ratio’s are not within the recommended range. Therefore, we will proceed with the default tresholds for the ligand activity analysis
After the ligand activity prediction, we will also infer the predicted target genes of these ligands in each contrast. For this ligand-target inference procedure, we also need to select which top n of the predicted target genes will be considered (here: top 250 targets per ligand). This parameter will not affect the ligand activity predictions. It will only affect ligand-target visualizations and construction of the intercellular regulatory network during the downstream analysis. We recommend users to test other settings in case they would be interested in exploring fewer, but more confident target genes, or vice versa.
top_n_target = 250
The NicheNet ligand activity analysis can be run in parallel for each receiver cell type, by changing the number of cores as defined here. Using more cores will speed up the analysis at the cost of needing more memory. This is only recommended if you have many receiver cell types of interest.
verbose = TRUE
cores_system = 8
n.cores = min(cores_system, celltype_de$cluster_id %>% unique() %>% length())
Running the ligand activity prediction will take some time (the more cell types and contrasts, the more time)
ligand_activities_targets_DEgenes = suppressMessages(suppressWarnings(
get_ligand_activities_targets_DEgenes(
receiver_de = celltype_de,
receivers_oi = intersect(receivers_oi, celltype_de$cluster_id %>% unique()),
ligand_target_matrix = ligand_target_matrix,
logFC_threshold = logFC_threshold,
p_val_threshold = p_val_threshold,
p_val_adj = p_val_adj,
top_n_target = top_n_target,
verbose = verbose,
n.cores = n.cores
)
))
You can check the output of the ligand activity and ligand-target inference here:
ligand_activities_targets_DEgenes$ligand_activities %>% head(20)
## # A tibble: 20 × 8
## # Groups: receiver, contrast [1]
## ligand activity contrast target ligand_target_weight receiver
## <chr> <dbl> <chr> <chr> <dbl> <chr>
## 1 A2M 0.0519 PreE-PreNE ACTA2 0.00715 CD4T
## 2 A2M 0.0519 PreE-PreNE ALOX5AP 0.00727 CD4T
## 3 A2M 0.0519 PreE-PreNE ARID5B 0.00719 CD4T
## 4 A2M 0.0519 PreE-PreNE BCL2L11 0.00797 CD4T
## 5 A2M 0.0519 PreE-PreNE BCL3 0.00825 CD4T
## 6 A2M 0.0519 PreE-PreNE BHLHE40 0.00945 CD4T
## 7 A2M 0.0519 PreE-PreNE BST2 0.00662 CD4T
## 8 A2M 0.0519 PreE-PreNE CDK6 0.00889 CD4T
## 9 A2M 0.0519 PreE-PreNE CDKN2A 0.00716 CD4T
## 10 A2M 0.0519 PreE-PreNE CDKN2C 0.00726 CD4T
## 11 A2M 0.0519 PreE-PreNE CKS1B 0.00744 CD4T
## 12 A2M 0.0519 PreE-PreNE CORO1C 0.00677 CD4T
## 13 A2M 0.0519 PreE-PreNE CSF2 0.00859 CD4T
## 14 A2M 0.0519 PreE-PreNE DDX60 0.00714 CD4T
## 15 A2M 0.0519 PreE-PreNE DUSP4 0.00677 CD4T
## 16 A2M 0.0519 PreE-PreNE DUSP5 0.00836 CD4T
## 17 A2M 0.0519 PreE-PreNE FKBP5 0.00723 CD4T
## 18 A2M 0.0519 PreE-PreNE GADD45G 0.00768 CD4T
## 19 A2M 0.0519 PreE-PreNE GAPDH 0.00729 CD4T
## 20 A2M 0.0519 PreE-PreNE GEM 0.00689 CD4T
## # ℹ 2 more variables: direction_regulation <fct>, activity_scaled <dbl>
In the previous steps, we calculated expression, differential expression and NicheNet ligand activity. In the final step, we will now combine all calculated information to rank all sender-ligand—receiver-receptor pairs according to group/condition specificity. We will use the following criteria to prioritize ligand-receptor interactions:
We will combine these prioritization criteria in a single aggregated prioritization score. In the default setting, we will weigh each of these criteria equally (scenario = "regular"). This setting is strongly recommended. However, we also provide some additional setting to accomodate different biological scenarios. The setting scenario = "lower_DE" halves the weight for DE criteria and doubles the weight for ligand activity. This is recommended in case your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). The setting scenario = "no_frac_LR_expr" ignores the criterion “Sufficiently high expression levels of ligand and receptor in many samples of the same group”. This may be interesting for users that have data with a limited number of samples and don’t want to penalize interactions if they are not sufficiently expressed in some samples.
Finally, we still need to make one choice. For NicheNet ligand activity we can choose to prioritize ligands that only induce upregulation of target genes (ligand_activity_down = FALSE) or can lead potentially lead to both up- and downregulation (ligand_activity_down = TRUE). The benefit of ligand_activity_down = FALSE is ease of interpretability: prioritized ligand-receptor pairs will be upregulated in the condition of interest, just like their target genes. ligand_activity_down = TRUE can be harder to interpret because target genes of some interactions may be upregulated in the other conditions compared to the condition of interest. This is harder to interpret, but may help to pick up interactions that can also repress gene expression.
Here we will choose for setting ligand_activity_down = FALSE and focus specifically on upregulating ligands. At the end of this tutorial, we will explore the effect of setting ligand_activity_down = TRUE.
ligand_activity_down = FALSE
sender_receiver_tbl = sender_receiver_de %>% distinct(sender, receiver)
metadata_combined = SummarizedExperiment::colData(sce) %>% tibble::as_tibble()
if(!is.na(batches)){
grouping_tbl = metadata_combined[,c(sample_id, group_id, batches)] %>%
tibble::as_tibble() %>% distinct()
colnames(grouping_tbl) = c("sample","group",batches)
} else {
grouping_tbl = metadata_combined[,c(sample_id, group_id)] %>%
tibble::as_tibble() %>% distinct()
colnames(grouping_tbl) = c("sample","group")
}
prioritization_tables = suppressMessages(generate_prioritization_tables(
sender_receiver_info = abundance_expression_info$sender_receiver_info,
sender_receiver_de = sender_receiver_de,
ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes,
contrast_tbl = contrast_tbl,
sender_receiver_tbl = sender_receiver_tbl,
grouping_tbl = grouping_tbl,
scenario = "regular", # all prioritization criteria will be weighted equally
fraction_cutoff = fraction_cutoff,
abundance_data_receiver = abundance_expression_info$abundance_data_receiver,
abundance_data_sender = abundance_expression_info$abundance_data_sender,
ligand_activity_down = ligand_activity_down
))
Check the output tables
First: group-based summary table
prioritization_tables$group_prioritization_tbl %>% head(20)
## # A tibble: 20 × 18
## contrast group sender receiver ligand receptor lr_interaction id
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 PreE-PreNE PreE CD4T macrophages BTLA TNFRSF14 BTLA_TNFRSF14 BTLA…
## 2 PreE-PreNE PreE Fibroblast macrophages CSF1 SIRPA CSF1_SIRPA CSF1…
## 3 PreE-PreNE PreE macrophages CD4T IL15 IL2RB IL15_IL2RB IL15…
## 4 PreE-PreNE PreE CD4T macrophages IFNG IFNGR2 IFNG_IFNGR2 IFNG…
## 5 PreE-PreNE PreE macrophages CD4T CD274 PDCD1 CD274_PDCD1 CD27…
## 6 PreE-PreNE PreE macrophages CD4T CD86 CTLA4 CD86_CTLA4 CD86…
## 7 PreE-PreNE PreE macrophages macrophages MMP9 CD44 MMP9_CD44 MMP9…
## 8 PreE-PreNE PreE macrophages macrophages IL15 IL2RG IL15_IL2RG IL15…
## 9 PreE-PreNE PreE macrophages CD4T NECTI… TIGIT NECTIN2_TIGIT NECT…
## 10 PreE-PreNE PreE macrophages CD4T IL15 IL2RG IL15_IL2RG IL15…
## 11 PreE-PreNE PreE Fibroblast macrophages RARRE… NRP2 RARRES1_NRP2 RARR…
## 12 PreE-PreNE PreE macrophages CD4T SIRPA TIGIT SIRPA_TIGIT SIRP…
## 13 PreE-PreNE PreE Fibroblast macrophages BST2 LILRA5 BST2_LILRA5 BST2…
## 14 PreE-PreNE PreE Fibroblast macrophages BST2 LILRB3 BST2_LILRB3 BST2…
## 15 PreE-PreNE PreE macrophages CD4T CD48 PDCD1 CD48_PDCD1 CD48…
## 16 PreE-PreNE PreE macrophages macrophages SIRPA CD47 SIRPA_CD47 SIRP…
## 17 PreE-PreNE PreE macrophages macrophages MMP9 IFNAR1 MMP9_IFNAR1 MMP9…
## 18 PreE-PreNE PreE macrophages macrophages LGALS3 ANXA2 LGALS3_ANXA2 LGAL…
## 19 PreE-PreNE PreE CD4T Fibroblast IFNG IFNGR2 IFNG_IFNGR2 IFNG…
## 20 PreE-PreNE PreE macrophages macrophages MMP9 ITGAM MMP9_ITGAM MMP9…
## # ℹ 10 more variables: scaled_lfc_ligand <dbl>,
## # scaled_p_val_ligand_adapted <dbl>, scaled_lfc_receptor <dbl>,
## # scaled_p_val_receptor_adapted <dbl>, max_scaled_activity <dbl>,
## # scaled_pb_ligand <dbl>, scaled_pb_receptor <dbl>,
## # fraction_expressing_ligand_receptor <dbl>, prioritization_score <dbl>,
## # top_group <chr>
This table gives the final prioritization score of each interaction, and the values of the individual prioritization criteria.
With this step, all required steps are finished. Now, we can optionally still run the following steps * Calculate the across-samples expression correlation between ligand-receptor pairs and target genes * Prioritize communication patterns involving condition-specific cell types through an alternative prioritization scheme
Here we will only focus on the expression correlation step:
In multi-sample datasets, we have the opportunity to look whether expression of ligand-receptor across all samples is correlated with the expression of their by NicheNet predicted target genes. This is what we will do with the following line of code:
lr_target_prior_cor = lr_target_prior_cor_inference(
receivers_oi = prioritization_tables$group_prioritization_tbl$receiver %>% unique(),
abundance_expression_info = abundance_expression_info,
celltype_de = celltype_de,
grouping_tbl = grouping_tbl,
prioritization_tables = prioritization_tables,
ligand_target_matrix = ligand_target_matrix,
logFC_threshold = logFC_threshold,
p_val_threshold = p_val_threshold,
p_val_adj = p_val_adj
)
To avoid needing to redo the analysis later, we will here to save an output object that contains all information to perform all downstream analyses.
path = "./"
multinichenet_output = list(
celltype_info = abundance_expression_info$celltype_info,
celltype_de = celltype_de,
sender_receiver_info = abundance_expression_info$sender_receiver_info,
sender_receiver_de = sender_receiver_de,
ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes,
prioritization_tables = prioritization_tables,
grouping_tbl = grouping_tbl,
lr_target_prior_cor = lr_target_prior_cor
)
multinichenet_output = make_lite_output(multinichenet_output)
save = FALSE
if(save == TRUE){
saveRDS(multinichenet_output, paste0(path, "multinichenet_output.rds"))
}
We suggest to split up the analysis in at least two scripts: the code to create this MultiNicheNet output object, and the code to analyze and interpret this output. For sake of demonstration, we will continue here in this vignette.
In a first instance, we will look at the broad overview of prioritized interactions via condition-specific Chordiagram circos plots. The aim of this visualizatin is to provide a summary of the top prioritized senderLigand-receiverReceptor interactions per condition (between all cell types or between cell type pairs of interest).
We will look here at the top 50 predictions across all contrasts, senders, and receivers of interest.
prioritized_tbl_oi_all = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
top_n = 50,
rank_per_group = FALSE
)
prioritized_tbl_oi =
multinichenet_output$prioritization_tables$group_prioritization_tbl %>%
filter(id %in% prioritized_tbl_oi_all$id) %>%
distinct(id, sender, receiver, ligand, receptor, group) %>%
left_join(prioritized_tbl_oi_all)
prioritized_tbl_oi$prioritization_score[is.na(prioritized_tbl_oi$prioritization_score)] = 0
senders_receivers = union(prioritized_tbl_oi$sender %>% unique(), prioritized_tbl_oi$receiver %>% unique()) %>% sort()
colors_sender = RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral') %>% magrittr::set_names(senders_receivers)
colors_receiver = RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral') %>% magrittr::set_names(senders_receivers)
circos_list = make_circos_group_comparison(prioritized_tbl_oi, colors_sender, colors_receiver)
Whereas these ChordDiagram circos plots show the most specific interactions per group, they don’t give insights into the data behind these predictions. Because inspecting the data behind the prioritization is recommended to decide on which interactions to validate, we created several functionalities to do this.
Therefore we will now generate “interpretable bubble plots” that indicate the different prioritization criteria used in MultiNicheNet.
In the next type of plots, we will visualize the following prioritization criteria used in MultiNicheNet: * 1) differential expression of ligand and receptor: the per-sample scaled product of normalized ligand and receptor pseudobulk expression * 2) the scaled ligand activities * 3) cell-type specificity of ligand and receptor.
As a further help for users to further prioritize, we also visualize: * the condition-average of the fraction of cells expressing the ligand and receptor in the cell types of interest * the level of curation of these LR pairs as defined by the Intercellular Communication part of the Omnipath database (https://omnipathdb.org/)
We will create this plot for BRCA group specific interactions of the overall top50 interactions that we visualized in the Circos Chorddiagrams above:
group_oi = "PreE"
prioritized_tbl_oi_group1_50 = prioritized_tbl_oi_all %>%
filter(group == group_oi)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi_group1_50 %>% inner_join(lr_network_all)
)
plot_oi
Some notes about this plot:
* Samples that were left out of the DE analysis (because too few cells in that celltype-sample combination) are indicated with a smaller dot. This helps to indicate the samples that did not contribute to the calculation of the logFC, and thus not contributed to the final prioritization.
* As you can see, interactions like the SIGLEC1-CD47 interaction do not have Omnipath DB scores. This is because this LR pair was not documented by the Omnipath LR database. Instead it was documented by the Verschueren database as can be seen in the table (
lr_network_all %>% filter(ligand == "SIGLEC1" & receptor == "CD47")).
Question: which interactions seem most relevant to you to explore for further validation?
–
We encourage users to make these plots also for the other groups, like we will do now first for the PreNE group
group_oi = "PreNE"
prioritized_tbl_oi_group2_50 = prioritized_tbl_oi_all %>%
filter(group == group_oi)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi_group2_50 %>% inner_join(lr_network_all)
)
plot_oi
Question: how could you interpret the difference in TGFB1 vs TGFB3 ligand activity here?
–
As you could observe from the Circos ChordDiagram and Interpretable Bubble plots above: we find more specific interactions for the PreE group than for the PreNE group here.
If you want to visualize more interactions specific for a group of interest, so not restricted to e.g. the top50 overall, but the top50 for a group of interest, you can run the following:
prioritized_tbl_oi_group2_50 = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
50,
groups_oi = group_oi)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi_group2_50 %>% inner_join(lr_network_all)
)
plot_oi
Typically, there are way more than 50 differentially expressed and active ligand-receptor pairs per group across all sender-receiver combinations. Therefore it might be useful to zoom in on specific cell types as senders/receivers:
We will illustrate this for the “CD4T” cell type as receiver in the PreE group:
group_oi = "PreE"
prioritized_tbl_oi_group1_50 = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
50,
groups_oi = group_oi,
receivers_oi = "CD4T"
)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi_group1_50 %>% inner_join(lr_network_all)
)
plot_oi
And now as sender:
prioritized_tbl_oi_group1_50 = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
50,
groups_oi = group_oi,
senders_oi = "CD4T")
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi_group1_50 %>% inner_join(lr_network_all))
plot_oi
These two types of plots created above (Circos ChordDiagram and Interpretable Bubble Plot) for the most strongly prioritized interactions are the types of plot you should always create and inspect as an end-user.
The plots that we will discuss in the rest of the vignette are more optional, and can help to dive more deeply in the data. They are however not as necessary as the plots above.
So, let’s now continue with more detailed plots and downstream functionalities:
In the plots above, we showed some of the prioritized interactions, and focused on their expression and activity. These interactions were visualized as independent interactions. However, they are likely not functioning independently in a complex multicellular biological system: cells can send signals to other cells, who as a response to these signals produce extracellular signals themselves to give feedback to the original sender cells, or to propogate the signal to other cell types (“cascade”). In other words: ligands from cell type A may induce the expression of ligands and receptors in cell type B. These ligands and receptors can then be involved in other interactions towards cell type A and interactions towards cell type C. Etc.
Because one of the elements of MultiNicheNet is the ligand activity and ligand-target inference part of NicheNet, we can actually infer the predicted ligand/receptor-encoding target genes of prioritized ligand-receptor interactions. And as a result, we can get this type of functional insight in the biological system of interest, which we will demonstrate now.
First, we will showcase how to do this by considering target genes supported by NicheNet’s prior knowledge solely
First: get the target genes of prioritized ligand-receptor pairs (here focused on the overall top50 prioritized LR pairs that were visualized in the Circos ChordDiagrams above based on the prioritized_tbl_oi_all data frame)
lr_target_prior = prioritized_tbl_oi_all %>% inner_join(
multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>%
distinct(ligand, target, direction_regulation, contrast) %>% inner_join(contrast_tbl) %>% ungroup()
)
lr_target_df = lr_target_prior %>% distinct(group, sender, receiver, ligand, receptor, id, target, direction_regulation)
Second, subset on ligands/receptors as target genes
lr_target_df %>% filter(target %in% union(lr_network$ligand, lr_network$receptor))
## # A tibble: 719 × 8
## group sender receiver ligand receptor id target direction_regulation
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <fct>
## 1 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… ADM up
## 2 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… CCL5 up
## 3 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… CD44 up
## 4 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… CSF1 up
## 5 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… CXCL10 up
## 6 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… FAS up
## 7 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… IL32 up
## 8 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… ITGAL up
## 9 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… MMP9 up
## 10 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TN… SDC4 up
## # ℹ 709 more rows
Whereas these code blocks are just to demonstrate that this type of information is available in MultiNicheNet, the next block of code will infer the systems-wide intercellular regulatory network automatically:
network = infer_intercellular_regulatory_network(lr_target_df, prioritized_tbl_oi_all)
network$links %>% head()
## # A tibble: 6 × 6
## sender_ligand receiver_target direction_regulation group type weight
## <chr> <chr> <fct> <chr> <chr> <dbl>
## 1 CD4T_BTLA macrophages_MMP9 up PreE Ligand-T… 1
## 2 Fibroblast_CSF1 macrophages_MMP9 up PreE Ligand-T… 1
## 3 macrophages_IL15 CD4T_CSF2 up PreE Ligand-T… 1
## 4 macrophages_IL15 CD4T_IFNG up PreE Ligand-T… 1
## 5 CD4T_IFNG macrophages_CD274 up PreE Ligand-T… 1
## 6 CD4T_IFNG macrophages_CXCL9 up PreE Ligand-T… 1
network$nodes %>% head()
## # A tibble: 6 × 4
## node celltype gene type_gene
## <chr> <chr> <chr> <chr>
## 1 CD4T_BTLA CD4T BTLA ligand/receptor
## 2 macrophages_SIRPA macrophages SIRPA ligand/receptor
## 3 macrophages_CD47 macrophages CD47 ligand/receptor
## 4 CD4T_SIRPG CD4T SIRPG ligand/receptor
## 5 macrophages_TNFRSF14 macrophages TNFRSF14 ligand/receptor
## 6 Fibroblast_CSF1 Fibroblast CSF1 ligand
And this network can be visualized here in R by running:
colors_sender["Fibroblast"] = "pink" # the original yellow background with white font is not very readable
network_graph = visualize_network(network, colors_sender)
network_graph$plot
As you can see here: we can see see here that several prioritized ligands seem to be regulated by other prioritized ligands! But, it may be challenging sometimes to discern individual links when several interactions are shown. Therefore, inspection of the underlying data tables (network$links and network$nodes) may be necessary to discern individual interactions. It is also suggested to export these data tables into more sophisticated network visualization tools (e.g., CytoScape) for better inspection of this network.
To inspect interactions involving specific ligands, such as IFNG as example, we can run the following code:
network$nodes %>% filter(gene == "IFNG")
## # A tibble: 1 × 4
## node celltype gene type_gene
## <chr> <chr> <chr> <chr>
## 1 CD4T_IFNG CD4T IFNG ligand
IFNG as regulating ligand:
network$links %>% filter(sender_ligand == "CD4T_IFNG" & direction_regulation == "up" & group == "PreE")
## # A tibble: 9 × 6
## sender_ligand receiver_target direction_regulation group type weight
## <chr> <chr> <fct> <chr> <chr> <dbl>
## 1 CD4T_IFNG macrophages_CD274 up PreE Ligand-T… 1
## 2 CD4T_IFNG macrophages_CXCL9 up PreE Ligand-T… 1
## 3 CD4T_IFNG macrophages_HLA.G up PreE Ligand-T… 1
## 4 CD4T_IFNG macrophages_IL15 up PreE Ligand-T… 1
## 5 CD4T_IFNG macrophages_PDCD1LG2 up PreE Ligand-T… 1
## 6 CD4T_IFNG Fibroblast_BST2 up PreE Ligand-T… 1
## 7 CD4T_IFNG Fibroblast_CSF1 up PreE Ligand-T… 1
## 8 CD4T_IFNG Fibroblast_RARRES1 up PreE Ligand-T… 1
## 9 CD4T_IFNG macrophages_CD47 up PreE Ligand-T… 1
IFNG as regulated target:
network$links %>% filter(receiver_target == "CD4T_IFNG" & direction_regulation == "up" & group == "PreE")
## # A tibble: 15 × 6
## sender_ligand receiver_target direction_regulation group type weight
## <chr> <chr> <fct> <chr> <chr> <dbl>
## 1 macrophages_IL15 CD4T_IFNG up PreE Ligan… 1
## 2 macrophages_CD274 CD4T_IFNG up PreE Ligan… 1
## 3 macrophages_CD86 CD4T_IFNG up PreE Ligan… 1
## 4 macrophages_NECTIN2 CD4T_IFNG up PreE Ligan… 1
## 5 macrophages_SIRPA CD4T_IFNG up PreE Ligan… 1
## 6 macrophages_CD48 CD4T_IFNG up PreE Ligan… 1
## 7 CD4T_BTLA CD4T_IFNG up PreE Ligan… 1
## 8 macrophages_CD47 CD4T_IFNG up PreE Ligan… 1
## 9 macrophages_TNFRSF14 CD4T_IFNG up PreE Ligan… 1
## 10 macrophages_B2M CD4T_IFNG up PreE Ligan… 1
## 11 macrophages_LGALS3 CD4T_IFNG up PreE Ligan… 1
## 12 macrophages_PDCD1LG2 CD4T_IFNG up PreE Ligan… 1
## 13 macrophages_CXCL9 CD4T_IFNG up PreE Ligan… 1
## 14 macrophages_LYZ CD4T_IFNG up PreE Ligan… 1
## 15 Fibroblast_VCAM1 CD4T_IFNG up PreE Ligan… 1
Ligand- and receptor-encoding target genes that were shown here are predicted as target genes of ligands based on prior knowledge. However, it is uncertain whether they are also potentially active in the system under study: e.g., it is possible that some genes are regulated by their upstream ligand only in cell types that are not studied in this context. To increase the chance that inferred ligand-target links are potentially active, we can use the multi-sample nature of this data to filter target genes based on expression correlation between the upstream ligand-receptor pair and the downstream target gene. This is under the assumption that target genes that show across-sample expression correlation with their upstream ligand-receptor pairs may be more likely to be true active target genes than target genes that don’t show this pattern. This correlation was calculated in the (optional) step 7 of the MultiNicheNet analysis.
In the next subsection of the inference of intercellular regulator networks, we will showcase how to consider target genes that are both supported by NicheNet’s prior knowledge and expression correlation.
Now, we will filter out correlated ligand-receptor –> target links that both show high expression correlation (pearson correlation > 0.33 in this example) and have some prior knowledge to support their link.
lr_target_prior_cor_filtered =
multinichenet_output$prioritization_tables$group_prioritization_tbl$group %>% unique() %>%
lapply(function(group_oi){
lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>%
inner_join(
multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>%
distinct(ligand, target, direction_regulation, contrast)
) %>%
inner_join(contrast_tbl) %>% filter(group == group_oi)
lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>%
filter(direction_regulation == "up") %>%
filter( (rank_of_target < top_n_target) & (pearson > 0.33))
lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>%
filter(direction_regulation == "down") %>%
filter( (rank_of_target < top_n_target) & (pearson < -0.33))
lr_target_prior_cor_filtered = bind_rows(
lr_target_prior_cor_filtered_up,
lr_target_prior_cor_filtered_down
)
}) %>% bind_rows()
lr_target_df = lr_target_prior_cor_filtered %>%
distinct(group, sender, receiver, ligand, receptor, id, target, direction_regulation)
network = infer_intercellular_regulatory_network(lr_target_df, prioritized_tbl_oi_all)
network$links %>% head()
## # A tibble: 6 × 6
## sender_ligand receiver_target direction_regulation group type weight
## <chr> <chr> <fct> <chr> <chr> <dbl>
## 1 macrophages_LGALS3 macrophages_MMP9 up PreE Ligand… 1
## 2 macrophages_MMP9 macrophages_MMP9 up PreE Ligand… 1
## 3 macrophages_ICAM1 macrophages_MMP9 up PreE Ligand… 1
## 4 macrophages_LGALS9 macrophages_MMP9 up PreE Ligand… 1
## 5 CD4T_IFNG macrophages_CD274 up PreE Ligand… 1
## 6 CD4T_IFNG macrophages_CXCL9 up PreE Ligand… 1
network$nodes %>% head()
## # A tibble: 6 × 4
## node celltype gene type_gene
## <chr> <chr> <chr> <chr>
## 1 macrophages_SIRPA macrophages SIRPA ligand/receptor
## 2 CD4T_SIRPG CD4T SIRPG ligand/receptor
## 3 CD4T_BTLA CD4T BTLA ligand/receptor
## 4 macrophages_CD47 macrophages CD47 ligand/receptor
## 5 macrophages_TNFRSF14 macrophages TNFRSF14 ligand/receptor
## 6 macrophages_LGALS3 macrophages LGALS3 ligand
network_graph = visualize_network(network, colors_sender)
network_graph$plot
As can be expected, we see fewer links here than in the previously generated intercellular regulatory network. The links that are not present anymore in this network are those ligand-target links that are not supported by high across-sample expression correlation. In conclusion, the links visualized here are the most trustworthy ones, since they are both supported by prior knowledge and expression correlation.
Interestingly, ligands/receptors visualized in this network can be considered as additionally prioritized because they are not only a prioritized ligand/receptor but also a target gene of another prioritized ligand-receptor interaction! So, we can also use this network to further prioritize differential CCC interactions. We can get these interactions as follows:
network$prioritized_lr_interactions
## # A tibble: 44 × 5
## group sender receiver ligand receptor
## <chr> <chr> <chr> <chr> <chr>
## 1 PreE macrophages macrophages LGALS3 ANXA2
## 2 PreE macrophages macrophages MMP9 ITGB2
## 3 PreE macrophages macrophages MMP9 CD44
## 4 PreE macrophages macrophages ICAM1 ITGB2
## 5 PreE macrophages macrophages MMP9 IFNAR1
## 6 PreE macrophages macrophages LGALS9 CD44
## 7 PreE CD4T macrophages IFNG IFNGR2
## 8 PreE macrophages macrophages MMP9 ITGAM
## 9 PreE macrophages macrophages LGALS9 CD47
## 10 PreE macrophages macrophages SIRPA CD47
## # ℹ 34 more rows
prioritized_tbl_oi_network = prioritized_tbl_oi_all %>% inner_join(
network$prioritized_lr_interactions)
prioritized_tbl_oi_network
## # A tibble: 44 × 8
## group sender receiver ligand receptor id prioritization_score
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 PreE CD4T macrophages BTLA TNFRSF14 BTLA_TNF… 0.949
## 2 PreE Fibroblast macrophages CSF1 SIRPA CSF1_SIR… 0.944
## 3 PreE macrophages CD4T IL15 IL2RB IL15_IL2… 0.937
## 4 PreE CD4T macrophages IFNG IFNGR2 IFNG_IFN… 0.934
## 5 PreE macrophages CD4T CD274 PDCD1 CD274_PD… 0.932
## 6 PreE macrophages CD4T CD86 CTLA4 CD86_CTL… 0.921
## 7 PreE macrophages macrophages MMP9 CD44 MMP9_CD4… 0.918
## 8 PreE macrophages CD4T NECTIN2 TIGIT NECTIN2_… 0.910
## 9 PreE macrophages CD4T IL15 IL2RG IL15_IL2… 0.907
## 10 PreE macrophages CD4T SIRPA TIGIT SIRPA_TI… 0.905
## # ℹ 34 more rows
## # ℹ 1 more variable: prioritization_rank <dbl>
Visualize now the expression and activity of these interactions for the PreE group
group_oi = "PreE"
prioritized_tbl_oi_group1 = prioritized_tbl_oi_network %>% filter(group == group_oi)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi_group1 %>% inner_join(lr_network_all)
)
plot_oi
To summarize: this interpretable bubble plot is an important and helpful plot because: 1) these LR interactions are all in the overall top50 of condition-specific interactions 2) they are a likely interaction inducing one or more other prioritized LR interaction and/or they are regulated by one or more other prioritized LR interactions. Because of this, interactions in this plot may be interesting candidates for follow-up experimental validation.
Note: These networks were generated by only looking at the top50 interactions overall. In practice, we encourage users to explore more hits than the top50, certainly if many cell type pairs are considered in the analysis.
All the previous were informative for interactions where both the sender and receiver cell types are captured in the data and where ligand and receptor are sufficiently expressed at the RNA level. However, these two conditions are not always fulfilled and some interesting cell-cell communication signals may be missed as a consequence. Can we still have an idea about these potentially missed interactions? Yes, we can.
In the next type of plot, we plot all the ligand activities (both scaled and absolute activities) of each receiver-condition combination. This can give us some insights in active signaling pathways across conditions. Note that we can thus show top ligands based on ligand activity - irrespective and agnostic of expression in sender. Benefits of this analysis are the possibility to infer the activity of ligands that are expressed by cell types that are not in your single-cell dataset or that are hard to pick up at the RNA level.
The following block of code will show how to visualize the activities for the top5 ligands for each receiver cell type - condition combination:
ligands_oi = multinichenet_output$prioritization_tables$ligand_activities_target_de_tbl %>%
inner_join(contrast_tbl) %>%
group_by(group, receiver) %>% filter(direction_regulation == "up") %>%
distinct(ligand, receiver, group, activity) %>%
top_n(5, activity) %>%
pull(ligand) %>% unique()
plot_oi = make_ligand_activity_plots(
multinichenet_output$prioritization_tables,
ligands_oi,
contrast_tbl,
widths = NULL)
plot_oi
Interestingly, we can here see a clear interferon signature among the upregulated genes in all cell types in the PreE patient group. The usefulness of this analysis: it can help you in having an idea about relevant ligands not captured in the data at hand but with a strong predicted target gene signature in one of the cell types in the data.
Note you can replace the automatically determined ligands_oi by any set of ligands that are of interest to you.
With this plot/downstream analysis, we end the overview of visualizations that can help you in finding interesting hypotheses about important differential ligand-receptor interactions in your data. In case you ended up with a shortlist of interactions for further checks and potential experimental validation, we recommend going over the visualizations that are introduced in the next section. They are some additional “sound checks” for your shortlist of interactions. However, we don’t recommend generating these plots before having thoroughly analyzed and inspected all the previous visualizations. Only go further now if you understood all the previous steps to avoid getting more overwhelmed.
Even though the interpretable bubble plots already provide a lot of information, they do not visualize the specific target genes downstream of the prioritized interactions. Hereby, we still miss some interesting functional information and we cannot assess whether high activity values may be due to a reasonable number of specific target genes or not. Therefore we will now go over some visualizations to inspect target genes downstream of prioritized ligand-receptor interactions.
In this type of plot, we can visualize the ligand activities for a group-receiver combination, and show the predicted ligand-target links, and also the expression of the predicted target genes across samples.
For this, we now need to define a receiver cell type of interest. As example, we will take CD4T cells as receiver, and look at the top 10 senderLigand-receiverReceptor pairs with these cells as receiver.
group_oi = "PreE"
receiver_oi = "CD4T"
prioritized_tbl_oi_group1_10 = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
10,
groups_oi = group_oi,
receivers_oi = receiver_oi)
combined_plot = make_ligand_activity_target_plot(
group_oi,
receiver_oi,
prioritized_tbl_oi_group1_10,
multinichenet_output$prioritization_tables,
multinichenet_output$ligand_activities_targets_DEgenes, contrast_tbl,
multinichenet_output$grouping_tbl,
multinichenet_output$celltype_info,
ligand_target_matrix,
plot_legend = FALSE)
combined_plot
## $combined_plot
##
## $legends
One observation we can make here is that several genes upregulated in the PreE-group are high-confident target genes of IL15 and CD274 (dark purple - high regulatory potential scores). Most of these genes are potential target genes of both these ligands, but some specific genes are present as well.
Whereas this plot just showed the top ligands for a certain receiver-contrast, you can also zoom in on specific ligands of interest. As example, we will look at “TGFB1”,“TGFB3”,“CTSD” (which came out as interesting ligand based on the sender-agnostic ligand activity procedure) in Fibroblasts in PreNE:
group_oi = "PreNE"
receiver_oi = "Fibroblast"
ligands_oi = c("TGFB1","TGFB3","CTSD")
prioritized_tbl_ligands_oi = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
100000,
groups_oi = group_oi,
receivers_oi = receiver_oi
) %>% filter(ligand %in% ligands_oi) # ligands should still be in the output tables of course
combined_plot = make_ligand_activity_target_plot(
group_oi,
receiver_oi,
prioritized_tbl_ligands_oi,
multinichenet_output$prioritization_tables,
multinichenet_output$ligand_activities_targets_DEgenes,
contrast_tbl,
multinichenet_output$grouping_tbl,
multinichenet_output$celltype_info,
ligand_target_matrix,
plot_legend = FALSE)
combined_plot
## $combined_plot
##
## $legends
In summary, these “Ligand activity - target gene combination plots” show how well ligand-target links are supported by general prior knowledge, but not whether they are likely to be active in the system under study. That’s what we will look at now.
In the previous plot, target genes were shown that are predicted as target gene of ligands based on prior knowledge. However, we can use the multi-sample nature of this data to filter target genes based on expression correlation between the upstream ligand-receptor pair and the downstream target gene. We will filter out correlated ligand-receptor –> target links that both show high expression correlation (pearson correlation > 0.33 in this example) and have some prior knowledge to support their link. Note that you can only make these visualization if you ran step 7 of the core MultiNicheNet analysis.
group_oi = "PreE"
receiver_oi = "CD4T"
lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>%
inner_join(
multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>%
distinct(ligand, target, direction_regulation, contrast)
) %>%
inner_join(contrast_tbl) %>% filter(group == group_oi, receiver == receiver_oi)
lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>%
filter(direction_regulation == "up") %>%
filter( (rank_of_target < top_n_target) & (pearson > 0.33)) # replace pearson by spearman if you want to filter on the spearman correlation
lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>%
filter(direction_regulation == "down") %>%
filter( (rank_of_target < top_n_target) & (pearson < -0.33)) # downregulation -- negative correlation - # replace pearson by spearman if you want to filter on the spearman correlation
lr_target_prior_cor_filtered = bind_rows(
lr_target_prior_cor_filtered_up,
lr_target_prior_cor_filtered_down)
Now we will visualize the top correlated target genes for the LR pairs that are also in the top 50 LR pairs discriminating the groups from each other:
prioritized_tbl_oi = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
50,
groups_oi = group_oi,
receivers_oi = receiver_oi)
lr_target_correlation_plot = make_lr_target_correlation_plot(
multinichenet_output$prioritization_tables,
prioritized_tbl_oi,
lr_target_prior_cor_filtered ,
multinichenet_output$grouping_tbl,
multinichenet_output$celltype_info,
receiver_oi,
plot_legend = FALSE)
lr_target_correlation_plot$combined_plot
This visualization can help users assess whether ligand-target links that are supported by general prior knowledge, are also potentially active in the system under study: target genes that show across-sample expression correlation with their upstream ligand-receptor pairs may be more likely true target genes than target genes that don’t show this pattern.
Even though this plot indicates the strength of the correlation between ligand-receptor expression and target gene expression, it’s hard to assess the pattern of correlation. To help users evaluate whether high correlation values are not due to artifacts, we provide the following LR-target expression scatter plot visualization for a selected LR pair and their targets:
ligand_oi = "CXCL16"
receptor_oi = "CXCR6"
sender_oi = "macrophages"
receiver_oi = "CD4T"
lr_target_scatter_plot = make_lr_target_scatter_plot(
multinichenet_output$prioritization_tables,
ligand_oi, receptor_oi, sender_oi, receiver_oi,
multinichenet_output$celltype_info,
multinichenet_output$grouping_tbl,
lr_target_prior_cor_filtered)
lr_target_scatter_plot
ligand_oi = "CD274"
receptor_oi = "PDCD1"
sender_oi = "macrophages"
receiver_oi = "CD4T"
lr_target_scatter_plot = make_lr_target_scatter_plot(
multinichenet_output$prioritization_tables,
ligand_oi, receptor_oi, sender_oi, receiver_oi,
multinichenet_output$celltype_info,
multinichenet_output$grouping_tbl,
lr_target_prior_cor_filtered)
lr_target_scatter_plot
The next type of “sound check” visualization will visualize potential signaling paths between ligands and target genes of interest. In addition to this visualization, we also get a network table documenting the underlying data source(s) behind each of the links shown in this graph. This analysis can help users to assess the trustworthiness of ligand-target predictions. This is strongly recommended before going into experimental validation of ligand-target links.
This inference of ‘prior knowledge’ ligand-receptor-to-target signaling paths is done similarly to the workflow described in the nichenetr package https://github.com/saeyslab/nichenetr/blob/master/vignettes/ligand_target_signaling_path.md
First read in the required networks:
if(organism == "human"){
sig_network = readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_human_21122021.rds")) %>%
mutate(from = make.names(from), to = make.names(to))
gr_network = readRDS(url("https://zenodo.org/record/7074291/files/gr_network_human_21122021.rds")) %>%
mutate(from = make.names(from), to = make.names(to))
ligand_tf_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_tf_matrix_nsga2r_final.rds"))
colnames(ligand_tf_matrix) = colnames(ligand_tf_matrix) %>% make.names()
rownames(ligand_tf_matrix) = rownames(ligand_tf_matrix) %>% make.names()
weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds"))
weighted_networks$lr_sig = weighted_networks$lr_sig %>% mutate(from = make.names(from), to = make.names(to))
weighted_networks$gr = weighted_networks$gr %>% mutate(from = make.names(from), to = make.names(to))
} else if(organism == "mouse"){
sig_network = readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_mouse_21122021.rds")) %>%
mutate(from = make.names(from), to = make.names(to))
gr_network = readRDS(url("https://zenodo.org/record/7074291/files/gr_network_mouse_21122021.rds")) %>%
mutate(from = make.names(from), to = make.names(to))
ligand_tf_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_tf_matrix_nsga2r_final_mouse.rds"))
colnames(ligand_tf_matrix) = colnames(ligand_tf_matrix) %>% make.names()
rownames(ligand_tf_matrix) = rownames(ligand_tf_matrix) %>% make.names()
weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final_mouse.rds"))
weighted_networks$lr_sig = weighted_networks$lr_sig %>% mutate(from = make.names(from), to = make.names(to))
weighted_networks$gr = weighted_networks$gr %>% mutate(from = make.names(from), to = make.names(to))
}
Define which ligand and target genes you want to focus on: let’s take a couple of CXCL16 target genes from the above scatter plots
ligand_oi = "CD274"
receptor_oi = "PDCD1"
targets_all = c("CCL4", "CXCL13","C16orf54")
active_signaling_network = nichenetr::get_ligand_signaling_path_with_receptor(
ligand_tf_matrix = ligand_tf_matrix,
ligands_all = ligand_oi,
receptors_all = receptor_oi,
targets_all = targets_all,
weighted_networks = weighted_networks,
top_n_regulators = 3
)
data_source_network = nichenetr::infer_supporting_datasources(
signaling_graph_list = active_signaling_network,
lr_network = lr_network %>% dplyr::rename(from = ligand, to = receptor),
sig_network = sig_network,
gr_network = gr_network
)
active_signaling_network_min_max = active_signaling_network
active_signaling_network_min_max$sig = active_signaling_network_min_max$sig %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
active_signaling_network_min_max$gr = active_signaling_network_min_max$gr %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
colors = c("ligand" = "purple", "receptor" = "orange", "target" = "royalblue", "mediator" = "grey60")
ggraph_signaling_path = make_ggraph_signaling_path(
active_signaling_network_min_max,
colors,
ligand_oi,
receptor_oi,
targets_all)
ggraph_signaling_path$plot
As mentioned, we can also inspect the network table documenting the underlying data source(s) behind each of the links shown in this graph. This analysis can help users to assess the trustworthiness of ligand-target predictions.
data_source_network %>% head()
## # A tibble: 6 × 5
## from to source database layer
## <chr> <chr> <chr> <chr> <chr>
## 1 AR CCL4 harmonizome_CHEA harmonizome_gr regulatory
## 2 AR CXCL13 HTRIDB HTRIDB regulatory
## 3 AR CXCL13 KnockTF KnockTF regulatory
## 4 CD274 C16orf54 CytoSig_all CytoSig regulatory
## 5 CD274 CCL4 CytoSig_all CytoSig regulatory
## 6 CD274 CXCL13 CytoSig_all CytoSig regulatory
The following type of “sound check” visualization will visualize the single-cell expression distribution of a ligand/receptor/target gene of interest in cell types of interest. This may be informative for users to inspect the data behind DE results. This can help users evaluate whether DE results at pseudobulk level were not due to artifacts.
Single-cell expression Violin plots of ligand-receptor interaction of interest: make_ligand_receptor_violin_plot
It is often useful to zoom in on specific ligand-receptor interactions of interest by looking in more detail to their expression at the single cell level
We will again check the CD274-PDCD1 interaction for sake of demonstration:
ligand_oi = "CD274"
receptor_oi = "PDCD1"
group_oi = "PreE"
sender_oi = "macrophages"
receiver_oi = "CD4T"
p_violin = make_ligand_receptor_violin_plot(
sce = sce,
ligand_oi = ligand_oi,
receptor_oi = receptor_oi,
group_oi = group_oi,
group_id = group_id,
sender_oi = sender_oi,
receiver_oi = receiver_oi,
sample_id = sample_id,
celltype_id = celltype_id)
p_violin
For the CD274 target genes we visualized the signaling paths for, we can also inspect their single-cell expression levels:
list_target_plots = lapply(targets_all, function(target_oi) {
p = make_target_violin_plot(sce = sce, target_oi = target_oi, receiver_oi = receiver_oi, group_oi = group_oi, group_id = group_id, sample_id, celltype_id = celltype_id)
})
list_target_plots
## [[1]]
##
## [[2]]
##
## [[3]]
Finally, we provide some visualizations to just inspect the DE results that were generated during the MultiNicheNet analysis.
group_oi = "PreE"
receiver_oi = "CD4T"
DE_genes = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>%
inner_join(contrast_tbl) %>%
filter(group == group_oi) %>%
arrange(p_val) %>%
filter(
receiver == receiver_oi &
logFC > 2 &
p_val <= 0.05 &
contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>%
pull(gene) %>% unique()
p_target = make_DEgene_dotplot_pseudobulk(
genes_oi = DE_genes,
celltype_info = multinichenet_output$celltype_info,
prioritization_tables = multinichenet_output$prioritization_tables,
celltype_oi = receiver_oi,
multinichenet_output$grouping_tbl)
p_target$pseudobulk_plot + ggtitle("DE genes (pseudobulk expression)")
p_target$singlecell_plot + ggtitle("DE genes (single-cell expression)")
Among these DE genes, you may be most interested in ligands or receptors
Ligands:
group_oi = "PreE"
receiver_oi = "CD4T"
DE_genes = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>%
inner_join(contrast_tbl) %>%
filter(group == group_oi) %>%
arrange(p_val) %>%
filter(
receiver == receiver_oi &
logFC > 1 &
p_val <= 0.05 &
contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>%
pull(gene) %>% unique()
DE_genes = DE_genes %>% intersect(lr_network$ligand)
p_target = make_DEgene_dotplot_pseudobulk(
genes_oi = DE_genes,
celltype_info = multinichenet_output$celltype_info,
prioritization_tables = multinichenet_output$prioritization_tables,
celltype_oi = receiver_oi,
multinichenet_output$grouping_tbl)
p_target$pseudobulk_plot + ggtitle("DE ligands (pseudobulk expression)")
p_target$singlecell_plot + ggtitle("DE ligands (single-cell expression)")
Receptors:
group_oi = "PreE"
receiver_oi = "CD4T"
DE_genes = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>%
inner_join(contrast_tbl) %>%
filter(group == group_oi) %>%
arrange(p_val) %>%
filter(
receiver == receiver_oi &
logFC > 1 &
p_val <= 0.05 &
contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>%
pull(gene) %>% unique()
DE_genes = DE_genes %>% intersect(lr_network$receptor)
p_target = make_DEgene_dotplot_pseudobulk(
genes_oi = DE_genes,
celltype_info = multinichenet_output$celltype_info,
prioritization_tables = multinichenet_output$prioritization_tables,
celltype_oi = receiver_oi,
multinichenet_output$grouping_tbl)
p_target$pseudobulk_plot + ggtitle("DE receptors (pseudobulk expression)")
p_target$singlecell_plot + ggtitle("DE receptors (single-cell expression)")
In the analysis, we chose for the setting ligand_activity_down = FALSE to focus specifically on upregulating ligands. Now, we will change this to ligand_activity_down = TRUE such that we explicitly prioritize ligands with high downregulatory activity as well.
ligand_activity_down = TRUE
prioritization_tables_alternative = suppressMessages(generate_prioritization_tables(
sender_receiver_info = abundance_expression_info$sender_receiver_info,
sender_receiver_de = sender_receiver_de,
ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes,
contrast_tbl = contrast_tbl,
sender_receiver_tbl = sender_receiver_tbl,
grouping_tbl = grouping_tbl,
scenario = "regular", # all prioritization criteria will be weighted equally
fraction_cutoff = fraction_cutoff,
abundance_data_receiver = abundance_expression_info$abundance_data_receiver,
abundance_data_sender = abundance_expression_info$abundance_data_sender,
ligand_activity_down = ligand_activity_down
))
prioritized_tbl_oi_all_alternative = get_top_n_lr_pairs(
prioritization_tables_alternative,
top_n = 50,
rank_per_group = FALSE
)
group_oi = "PreE"
prioritized_tbl_oi_group1_50_alternative = prioritized_tbl_oi_all_alternative %>%
filter(group == group_oi)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
prioritization_tables_alternative,
prioritized_tbl_oi_group1_50_alternative %>% inner_join(lr_network_all)
)
plot_oi
Now compare the unique interactions of this analysis versus the previous;
get top50 interactions of previous analysis with up-only:
prioritized_tbl_oi_all = get_top_n_lr_pairs(
multinichenet_output$prioritization_tables,
top_n = 50,
rank_per_group = FALSE
)
prioritized_tbl_oi_group1_50 = prioritized_tbl_oi_all %>%
filter(group == group_oi)
unique_up_down = prioritized_tbl_oi_group1_50_alternative$id %>% setdiff(prioritized_tbl_oi_group1_50$id)
unique_up = prioritized_tbl_oi_group1_50$id %>% setdiff(prioritized_tbl_oi_group1_50_alternative$id)
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
prioritization_tables_alternative,
prioritized_tbl_oi_group1_50_alternative %>% inner_join(lr_network_all) %>% filter(id %in% unique_up_down)
)
plot_oi
As expected, all of these interactions have a high value for scaled downregulatory activity. But they may be harder to interpret: for example, what does it mean that SEMA3C-PLXNA1 has a high downregulatory activity in PreE (and thus high upregulatory activity in PreNE), when the expression is not very clearly different between the two patient groups?
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
prioritization_tables,
prioritized_tbl_oi_group1_50 %>% inner_join(lr_network_all) %>% filter(id %in% unique_up)
)
plot_oi
plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
prioritization_tables,
prioritized_tbl_oi_group1_50 %>% inner_join(lr_network_all) %>% filter(!id %in% union(unique_up, unique_up_down))
)
plot_oi
Take-home message here:
The top interactions of both analyses settings seem very relevant and strongly differential. Some interactions are just ranked slightly higher/lower in one setting versus the other.
Choose the settings that are most appropriate to your dataset - choose the settings that will rank interactions according to how you will want to validate interactions.
This vignette covered all the details of a MultiNicheNet analysis, from running the algorithm to interpreting its output till the finest details. It’s important to realize that interpreting the output requires quite some time and different levels of iterations: start with the big picture and focus on the most differential LR pairs first. Then later, zoom in on target genes and perform the necessary “sound checks” when going further.